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FINAL  SCIENTIFIC  REPORT 
ONR  Contract  No.  N000U-76-C-0I82 

Introduction 

This  final  report  summarizes  our  research  activities  on  ONR 
Contract  No.  N00014-76-C-0182  from  February,  1976  to  January,  1983. 

''These  research  activities  resulted  in  twenty  four  publications 
and  twenty  two  oral  preserTEationsjiwere  in  part  supported  by  the 

NASA  and  the  AFOSR.  _ 

Fictitious  Gas  Design  Method 

The  design  technique  of  using  fictitious  gas  to  achieve  shock 
free  transonic  flight  conditions  had  been  studied  and  extended  to 
the  design  of  wings,  airfoils,  and  supersonic  conical  bodies.  Its 
practicality  and  the  ability  of  obtaining  near  optimal  designs  have 
been  summarized  in  the  review  paper  presented  at  the  International 
Congress  of  Aeronautical  Sciences  Meeting,  Seattle,  Washington. 

The  new  application  to  the  design  of  supersonic  conical  bodies 
was  reported  in  Dr.  Sritharan's  dissertation  and  also  mentioned  in 
the  ICAS  paper.  Only  the  example  of  a  circular  cone  at  an  angle 
of  incident  was  given.  Hopefully,  Dr.  Sritharan  will  be  able  to 
complete  this  study  at  ICAS  where  he  is  a  post-doctoral  fellow.  A 
general  conclusion  of  this  study  at  this  stage  is  that  the  applicability 
of  the  fictitious  method  in  conical  flows  is  more  limited  than  in 
plane  transonic  flow. 

Conical  Flows 

A  fully  conservative  finite  area  code  based  on  the  potential 


-2- 


approxlmatlon  for  conical  flows  has  been  developed  in  Dr.  Sritharan's 
dissertation.  This  work,  being  the  first  using  finite  area  approach 
and  one  that  generalizes  Jameson's  iteration  procedure  to  general 
curvilinear  systems,  was  presented  at  the  AIAA  Fluid,  Plasmadynamics 
and  Heat  Transfer  Conference,  St.  Louis,  Missouri  June,  1982.  An 
application  of  this  code  is  to  design  conical  bodies  at  incident 
with  shock  free  cross  flow,  mentioned  above. 

Unsteady  Transonic  Flow  Computations  with  Input  Pressure  Distributions 
Derived  from  our  earlier  belief  that  an  accurate  steady  pressure 
distribution  with  correct  shock  strength  and  location  is  important 
for  the  prediction  of  unsteady  transonic  responses,  an  inverse 
algorithm  IAF2,  accepting  pressure  distributions  instead  of  airfoil 
coordinates,  js  developed  to  provide  the  steady  flow  fields  needed 
in  unsteady  *,s/inputations.  This  algorithm  has  been  integrated  into 
our  version  of  LTRAN  and  the  time  linearized  code  UTFC  developed  in 
our  earlier  studies.  This  work  was  presented  at  the  AIAA  Fluid, 
Plasmadynamics  and  Heat  Transfer  Conference.  Very  good  agreement 
compared  with  experimental  results  in  predicting  phase  lags  and 
magnitudes  of  unsteady  responses  for  attached  transonic  flows  is 
obtained  using  this  procedure. 

Unsteady  Hind  Tunnel  Interference 

Alerted  by  previous  research  on  the  proper  far-field  boundary 
conditions  for  unsteady  transonic  flow,  we  have  concluded  in  a 
study,  as  part  of  Mr.  Przybytkowski's  dissertation,  that  two  dimensional 
wall  interferences  are  most  critical  in  the  low  to  moderate  reduced 
frequency  range;  that  resonance  conditions  can  be  predicted  using  linear 
acoustic  theory,  with  slight  modifications  due  to  nonlinear  effects, 
and  that  part  of  the  discrepancies  between  experimental  measurement 
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and  numerical  calculations  nw^y  be  attributed  to  the  uncertainities 
in  wall  conditions.  A  review  of  thi*  work  will  appear  in  a  book 
entitled,  "Recent  Advances  in  Numerical  Methods  in  Fluids,  Vol. 

IV,"  edited  by  Habashi.  A  copy  of  this  review  paper  is  in  Appendix 
C. 

Inviscid  Flows  with  Shock  Induced  Vorticity 

In  order  to  compute  the  flow  after  a  moderately  strong  shock 
accurately,  we  have  studied  and  confirmed  the  logartithmic  nature 
of  the  local  solution  near  the  shock  root  and  found  that  the  vorticity 
induced  by  the  curvature  of  the  shock,  despite  the  well  known 
singularity  at  the  root,  is  finite  and  given  by  a  formula  in  terms 
of  the  upstream  Mach  number,  density  and  the  curvature  of  the  body 
at  the  shock  root.  A  technical  note  on  this  result  is  published  in 
the  AIAA  Journal  attached  as  Appendix  D. 

Conclusion 

During  the  course  of  this  research,  we  have  completed  various 
studies  in  transonic  flow,  resulting  in  publications,  three  doctoral 
dissertations  and  two  master  reports,  and  seven  presentations.  These 
results  reflect  not  only  this  support,  but  also  AFOSR  support.  We 
conclude  this  report  with  a  list  of  publications  from  January,  1977 
through  June,  1983. 
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A  SIMPLE,  ACCURATE  AND  EFFICIENT  ALGORITHM  ■•OR  UNSTEADY 
TRANSONIC  FLOW 

K.-Y.  Fung 

Aerospace  and  Mechanical  Engineering 
University  of  Arizona 
Tucson,  Arizona 

1.  LNTRODUCTION 

This  chapter  describes  an  algorithm  for  computing  an 
airfoil's  response  to  small  unsteady  perturbations.  The 
development  of  the  algorithm  is  sufficiently  detailed  for 
the  experienced  'reader  to  code.  The  algorithm  is  simple, 
efficient,  and  incorporates  several  special  features  that 
make  it  an  effective  method  of  determining  the  aerodynamic 
response  of  airfoils  and,  using  strip  theory,  wings  with 
modest  sweepback.  The  special  features  include  an  accurate 
modeling  of  the  far-field  boundary  condition,  the  ability  to 
prescribe  the  airfoil's  steady-state  pressure  distribution 
and  determine  the  small  perturbation  steady-state  flow  field 
consistent  with  this  prescription,  and  the  capability  of 
genei^ting  results  for  low  reduced  frequencies  through  the 
superposition  of  the  results  for  a  single  step  change,  that 
is,  an  indicia  I  motion. 

We  begin  with  a  brief  discussion  of  the  importance  of 
these  flows,  then  state  the  governing  equations  and  discuss 
their  time  linearization.  Special  care  is  taken  to  describe 
the  treatment  of  moving  shock  waves.  The  equations  and 
boundary  conditions  are  then  discretized  using  the  usual  ADI 
scheme.  The  far-field  for  a  potential  vortex  with  wake  is 
used  to  provide  the  appropriate  far-field  boundary  condition. 
The  procedure  for  finding  Che  airfoil  chat  corresponds  to  a 
given  prescribed  pressure  distribution  is  then  described. 
The  advantages  of  using  an  Indiclal  response  are  then 
discussed.  Some  numerical  examples  and  a  discussion  of  the 
effects  of  wind  tunnel  walls  on  unsteady  transonic  testing 
follow.  Concluding  remarks  summarize  Che  capabllties  cf  the 
algorlchm  described. 
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2.  GOVERNING  EQUATIONS 

In  most  aeroelastic  problems,  strucCuril  deformations 
are  assumed  to  be  small  and,  hopefully,  will  remain  small 
through  aerodynamic  and  structural  damping.  However, 
oscillatory  modes  of  a  wing  and  its  control  surfaces  may 
extract  energy  from  free  stream  and  grow  until  the 
structure  disintegrates.  Flight  boundaries  that  prevent 
aircraft  from  operating  at  or  near  such  conditions  are 
established  in  aircraft  flight  tests.  An  essential  element 
in  aircraft  design  is  the  prediction  of  such  flutter 

boundaries.  To  establish  these  boundaries,  it  is  fair  to 
assume  that  the  unsteady  aerodynamic  disturbances  to  the 

already  established  steady  flow  field  about  the  aircraft  are 
small.  This  is  generally  true,  except  in  transonic  flow 
where  the  unsteady  pressure  fluctuations  may  be  large  due 
to  shock,  motions.  Nevertheless,  it  is  possible  to  correctly 
account  for  this  shoi.'  motion  by  linearizing  about  the 

steady-state  flow.  Care  must  be  taken  to  accurately  capture 
any  shock  waves  in  the  steady-state  flow. 

We  confine  our  study  to  thin  airfoils  with  a  thickness- 
to-chord  ratio  t  of  the  order  (1-Mi)^''"^.  Then  the  governing 
equation  for  the  perturbation  velocity  potential  -p  in  two 

space  dimensions  is 

-k'^M^i;itc~2kM^ipxfH  Y+l)M,|(?xl'?xx’^>yy  “  (D 

where  M*  is  the  free  stream  Mach  number,  y  is  the  ratio  of 
specific  heats,  and  k  is  the  reduced  frequency  obtained  from 
the  circular  frequency  nondimensionalized  by  the  free  stream 
velocity  U  divided  by  the  airfoil  chord  c.  This  equation, 

first*  derived  by  Lin  et  al.  [1],  can  be  obtained  from  a 
systematic  expansion  of  the  continuity  and  the  Bernoulli 

equations  in  the  thickness  ratio  t.  It  is  valid  for  inviscid 
and  Irrotatlonal  flow,  which  implies  an  attached  boundary 
layer,  and  the  potential  approximation  requires  that  any 

shock  wave  be  weak.  For  k  =  0(1)  the  nonlinear  term  is 

inconsequential  and  may  be  dropped. 

The  boundary  condition  on  the  airfoil  is  the  linearized 
flow  tangency  condition,  i.e., 

<|>y(x,0-,c)  -  t(yo  +  ^  (Y5{  +  kYg)|.  0  <  X  <  1,  (2) 

where  there  are  separate  functions  Y(x,  t)  for  the  upper  and 
lower  surface.  Here  the  airfoil  shape  has  been  decomposed 
into  a  fixed  part,  tY°, and  a  moving  part, eY“(x,t),  of  amplitude 
e.  In  the  wake  behind  the  trailing  edge  ,  the  pressure 
coef flclent ,Cp  ,1s  continuous  and  we  require 


j[  Ij  »  0,  across  y  =  0  for  x  >  1;  t,  i) 

here  denotes  jump  of  the  brackeCeij  quantities.  At 

large  distances,  all  steady  disturbances  should  decay 
properly  and  all  unsteady  disturbances  allowed  to  radiate  to 
infinity.  We  shall  discuss  the  implementation  of  these 
conditions  later. 

2.1  Far-Field  Boundary  Conditions 

At  large  distances  all  disturbances  must  decay. 
However,  the  size  of  the  coraputationa 1  domain  necessary  for 
trivial  boundary  conditions  has  a  profound  effect  on  t!ie 
efficiency  of  an  algorithm.  A  typical  icoustic  wave  of 
frequency  w  has  a  wave  length  a/,.-  or  ic/kM„),  and  thus 
requires  a  mesh  no  bigger  than  the  chorti  c  to  resolve  tnis 
wave  for  the  range  of  k  considered  here.  in  unsteady 
transonic  flow,  the  characteristic  time  scale  for  an 
unsteady  process  to  reach  a  new  equilibrium  point  is 
typically  the  time  required  for  the  airfoil  to  travel 
several  hundred  chord  lengths.  A  satisfactory  grid  needs  to 
be  100  chord  lengths  in  size  to  prevent  outgoing  waves  from 
reflecting  oft  the  grid  and  requires  about  130  intervals 
with  approximately  50  of  these  clustered  near  the  airfoil  to 
resolve  the  flow  details.  The  large  size  Insures  that  waves 
reflected  from  any  numerical  boundary  will  not  have 
sufficient  time  to  return  to  the  airfoil  and  contaminate  the 
results. 


Far  away,  a  change  of  incidence  of  the  airfoil  behaves 
like  a  potential  vortex  that  introduces  a  jump  in  potential, 
^r,  at  the  origin  at  time  tg.  Later,  the  disturbances  caused 
by  t|^s  change  propagate  in  time  and  in  space,  following  the 
characteristic  t  +  x  -  (x^  +  An  exact  solution  of 

the  linear  perturbation  potential  can  be  found  using 
analytical  techniques  [2].  In  two  dimensions  this  result  is 


Ar(to)  ^  g2 

”  — 2n~  f(x,0y,  (C-t„)) 


kM: 


where 


f(x,y,t)  -  H(t-hc-(x‘^+y^)^''-)  • 

f  ((t^+Zxt-y^)^^^  +  tj  l(t^+2xt-y^)^''‘ 

<  tan"^  -  +  tan“‘  - 

I  y  y 

a  -  (1-m2)1/2,  and, 

(0  t'  <  0 
H(t’)  -  { 

(1  t'  >  0 
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A  similar  result  can  be  found  in>  three  dimensions. 


For  arbitrary  changes  in  circulation,  r(c),  a 
superposition  of  this  result  gives  the  general  two- 

dimensional  far  field,  namely, 

1  3"  drcto) 

^>(x,y,t)  -  j  J  £(x,3y,  (t-to))  -  dt^  (^) 

o  ° 

For  a  given  circulation  history  Pit),  which  must  be  kept. 

Equation  (4)  can  be  evaluated  very  efficiently  for  each  point 
of  the  far-field  boundary.  Using  the  result  (4)  .  a 

satisfactory  computational  domain  can  be  as  small  as  15 

chords,  providing  nearly  a  factor  of  10  savings  in  the 

storage  requirements  and  execution  time. 

2.2  Time-Linearization  of  the  Governing  Equations 

An  effective  way  to  solve  Equation  il)  is  to  first 
obtain  the  steady-state  solution,  say  ■^^(x,  y),  satisfying  the 
steady  boundary  condition, 

0)  =.  tYO,  0  X  <  1, 

and  then  compute  which  will  deviate  only  slightly  from  po 
for  small  unsteady  disturbances,  i.e.,  1 4)  -  4>o  I  »  0(e),  except 
in  the  regions  of  shock  motion.  This  suggests  the 

introduction  of  an  unsteady  potential  v(x,  y,  t)  defined  by 

<t>(x,  y,  t)  -  4iO(x,  y)  +  et(x,  y,  t). 

Substituting  this  identity  into  Equation  (1)  and  neglecting 
terms  0( e^),  one  obtains  the  time-linearized  equation  for 
unsteady  small  disturbances: 

-k^M^Vtt  -  2kM^ii/xt  +  (ll-M^  -  (y  +  Vx}x  +  +yy  “  0. 

(5) 

The  corresponding  boundary  conditions  are 

♦yCx.  0,  t)  -  Y5J  +  kYjJ.  q<  X  <1;  (ba) 

[[I'x  +  k'l'cd  -  0,  I  <  X  (6b) 

It  should  be  noted  that  in  cases  where  the  shock  excursions 
are  large,  this  linearization  is  no  longer  possible.  The 

coefficient,  ( 1-M«  -  (Y+l)M»(>§l,in  Equation  (5)  obtained  from 

Che  steady  state  determines  where  the  flow  is  subsonic  or 
supersonic.  As  a  shock  moves,  the  region  ahead  must  be 
supersonic  and  the  region  behind,  subsonic  in  our 
approximation,  and  this  implies  that  the  spatial  properties 
must  change  accordingly.  Unless  the  region  where  the  shock 
travels  is  small.  Equation  (5)  is  inadequate  and  Equation  (1) 
should  be  used  instead.^ 


However,  using  Equation  (5),  once  the  steady  state  is 
determined  from  computations  or  experimental  measurnment s, 
the  geometry  of  the  airfoil  as  well  as  can  be  ignored 
completely  and  grid  refinements  usually  associ.ited  with  i 
complex  geometry  are  no  longer  needed. 

Equations  Cl)  and  (5)  can  be  further  simplified  by 
dropping  the  k2  terms;  as  we  noted  earlier  for  k  »  d(i),  the 
nonlinear  term  is  inconsequential.  Thus,  we  are  primarily 
interested  in  low-to-moderate  frequency  motions  ind  these 
are  the  most  likely  to  lead  to  flutter.  Indeed,  the 
appropriate  singular  limit  is  k  »  Thus  :he 

appropriate  equations  for  low  frequencies  are 

-2kMi^^^  +  [1-Mi  -  (Y-rl)Mi;>xl^yy  -^^vy  ’ 
in  place  of  Equation  (1),  and 

-2kMii;/xt  -  (y+1  Vx,x  »xy  '.’b) 

in  place  of  Equation  (5). 

2.3  Unsteady  Shock  Jump  Relations 

Equation  (7a)  admits  discontinuous  solutions,  that  is 
shocks,  satisfying  the  condition  derived  from  its 
conservative  form,  namely, 

dx 

-2kMi|^^]2^  -  [l-Mi-(Y  +  l)Mi5^l[;,^]^+||yi)^  .  J.  (8) 

together  with  the  condition  derived  from  irrotatlonal 
assum*ption 


Equation  (8)  states  that  the  shock  excursion, Xg(y,  t).can  be 
integrated  forward  in  time  knowing  the  jump  in  and  p.,. 
and  the  mean  ’■  'te,  across  the  discontinuity. 

Equation  (9)  descr  -e  shape  of  this  discontinuity  it  any 

instant  in  time. 

^or  thin  air  ransonic  Mach  numbers,  the  shock 

waves  that  occur  .  irly  normal  to  the  free  scream, 

therefore  the  jump  in  Vy  In  Equation  (8)  can  be  neglected, 
resulting  in  the  simplified  equation  for  Che  motion  of 
normal  shocks. 


<^xs  rH  / 
dt  “  2k  \ 


Mi-1  .  ) 
(Y+l)Mi 


Note  however,  Chat  Che  speed  of  the  shock  depends  on  and 
will  vary  with  y.  To  be  consistent  with  the  normal  shock 
assumption,  shock  motion  must  be  uniform  in  y  at  all  times. 
Thus,  Equation  (8)  is  only  integrated  at  one  y  location,  y  = 
y°,  for  each  shock  XgCc).  The  validity  of  this  approximation 
stems  from  the  observation  Chat  Che  shock  curvature  away 
from  the  shock  foot  where  it  is  singular  is  the  same  order 
as  Che  airfoil  curvature  and  hence  0(t/c)  and  consequently 
small.  A  careful  comparison  of  the  results  found  involving 
this  approximation  with  results  obtained  without  it,  (Yu  et 
al.  [3]),  confirmed  this  observation,  as  did  Fung  and  Chung 
(Al. 

Because  our  changes  Co  Che  steady  flow  field  are  0(e), 
it  is  clear  that  the  shock  motion  must  remain  small.  If  we 
linearize  the  shock  position  about  its  steady-state  value  x®. 


Xg  (t)  -  xg  +  5x(t), 


Chen  we  find 


and  we  note  that  for  the  reduced  frequencies  of  principal 
concern,  namely  k  •  we  must  have 

5  ^  (  ^273)  “  0(1), 

which  defines  Che  amplitudes  of  Che  unsteady  motion  for 
which  time  linearization  is  valid.  Thus,  for  k  ■  and 

with  ■  o(l),  we  may  use  the  time-linearized  result 

(7b). 


Near  Che  region  where  shock  motions  occur, 
perturbations  are  no  longer  small.  The  derivatives  of  if  can 
change  sign  abruptly  depending  on  whether  Che  shock  is  on 
Che  upstream  or  downstream  side.  This  would  invalidate  our 
small  disturbance  assumption  unless  if  is  created  differently 
to  account  for  large  disturbances  due  to  shock  motions. 
Consider  derivatives  of  the  potential  function  ^  at  time  t  - 
0.  The  potential  ^  will  have  a  discontinuous  slope  at  xg  as 
shown  in  Figure  1. 


Figure  1.  A  sketch  of  4°  and  f  in  the  vicinity  of  the 
unperturbed  shock. 

Let  us  assume  that  this  discontinuity  moves  to  a  new 
location  Xg  at  a  late  time  t  due  to  downstream  disturbances. 
To  account  for  this  change,  let  <p  be  expanded  about  the 

point  x|: 

•>(x,t)  ■  (f(xg,t)  +  (pxCxg.tXx-xl)  +  0(x-x§)2. 

If  <ii(x§,c)  is  replaced  by  v°(x|)  +  ev(x^,t),  and  we  require 
chat  <t)  be  continuous  at  Che  new  shock  point  Xg,  then  there 
must  be  no  Jump  in  at  that  point.  Thus, 

-  0  -  +  e{iiKx§,t)i|  +  {i<)>x°(xg)i](xg-xg) 

+  0(e(xs-x|))  +  0(xs-x|)2  +  .... 

Now,  since  both  Xg-x|  and  c  are  small,  the  change  of 

potential  at  Che  old  location  due  to  shock  motion  is  to  the 
lowest  order, 

--|<'S3(X3-xg), 

which  is  proportional  to  the  Jump  of  at  the  old  location 

and  Co  the  shock  excursion  (x^-xg).  To  accomplish  this 

shock  motion ,  we  analytically  (l.e.,  linearly)  continue 
upstream  to  the  steady'SCate  shock  position  where  we  have 
introduced  a  discontinuity  in  Thus  we  ignore  the  actual 
variation  in  in  regions  of  shock  motion.  This  allows  us  to 
account  for  shock  motions  provided  they  remain  small,  which 
will  generally  be  the  case  if  •  o(l). 


-  —  "er » 


3.  DlSCRtTIZATlUNS  OF  THE  GOVERNING  EQUATIONS 


a  » 


2kMi 

AtAxj  ’ 

2(l-ei) 

Ai-l/2,J 

Axi+i+Axi 

AXi 

2(l-ei) 

^1+1 / 2,  j 

Axi+i+Axi 

^xi+l 

2Ei-l 

Ai-l/2,j 

Axj+i+Axi 

Axi 

2^-1 

V2/3,j 

Axi+i+Axi 

Axi_i 

and 

Axi  -  Xi  -  Xi_l. 

Here  the  type-dependent  upwind  difference  scheme  of 
Murman  and  Cole  is  adopted  to  discretize  the  term  (AVx)x 
through  the  switching  function  This  scheme  reduces  to  a 

central-difference  scheme  if  the  flow  is  subsi>nic,  A^j  >  V, 
and  to  a  baqjkward-dif  ference  scheme  if  the  flow  is 
supersonic,  A^j  <  0;  in  both  cases,  the  matrix  Equation  (12) 
remains  diagonally  dominant. 

Across  the  shock  where  A^j  Is  positive  but  A^.j  j  is 
negative,  the  Jump  condition's  tinpieraented  by  setting  J  and 
e  to  zero  and  adding  (a+b)j['f>^jj  to  the  right-hand  side  of 
Equation  (12);  in  this  calculatioa 

and  (13) 

,n  ,n  r.n  1  ,n  n 

^n  ,  1  ,'t'lUs  :  n-I.js  -  lln,jsj  ^  n-l,js  -  Vjjs 

*  2  Axj^  ^^i-l 


Here  js  represents  the  location  y°  where  the  shock  speed  is 
computed,  normally  taken  to  be  zero. 


With  the  values  of  v"*"  determined,  the  new  values  of  v 
at  the  subsequent  time  level  calculated  using 


2kMi 

At 


'J'ij) 


or  the  algebraic  equatloa 


(-a,i+b+a,-bj 


’vij-il 


n+1 


Vij 

bij+lj 


Here 


civij  ~  V  “  2 


<Ayi+l  +  Ayi)Ayi’ 


(Ayi+i  +  Ayi)Ayi+i ’ 

2kMi 
AcAx^  ’ 


SyyVlj 


ia,-(a+t5),tjj 


Hj-i 

yij+ij 


and 

Ayi  -  yi  -  yi-1. 


(14) 


Again,  across  the  shock  location  the  jump  of  i>,  or  the 
term 


c 


^x(xi,y°) 


(15) 


is  added  to  the  right-hand  side  of  Equation  (14). 


The  full  procedure  is,  effectively. 


2kMi 

At 


T  ,  n+l  n  . 

^xUi,j  -  nj> 

°x^ij  +  J  '5yy('<'"lV+  '('ij) 


(16) 


with  the  jump  condition  Implemented  through  the  form  of 
Equations  (13)  and  (15).  It  can  be  shown  through  a  von 
Neumann  analysis  that,  shock  motion  aside,  this  scheme  is 
unconditionally  stable,  that  is,  there  is  no  restriction  on 
the  time  step  At.  However,  for  mixed  flow,  the  integration 
of  Equation  (11),  or  its  equivalent  Equation  (13),  for  the 
shock  motion  restricts  the  use  of  large  time  steps  for  both 

accuracy  and  stability.  _  _  ... _ 


I- 


3.1.  Boundary  Conditions 


A  typical  computational  domain  is  shown  in  Fig.  2.  On 
the  slit  representing  the  airfoil,  the  boundary  condition  (6a) 
becomes 


+  kY^.  0  <  XI  <  I 


and  in  the  wake  Equation  (6)  is  appro:<imated  by  using 
characteristics,  i.e.. 


which  is  unconditionally  stable  and  superior  to  the  implicit 
central-difference  scheme  or  the  upwind-backward  difference 
scheme  often  employed. 


Figure  2.  A  typical  grid  for  the  coaputatlons  described 
io  this  chapter. 


4.  PRESSURE  INPUT  STEADY  STATE 

The  solution  to  Equation  (7b)  depends  critically  on  the 
steady  state  that  is,  on  the  local  perturbed  pressure 
field.  This  is  especially  true  when  there  is  an  embedded 
supersonic  zone  where  exceeds  its  sonic  value 

(1'‘MJ)/(  In  this  case,  an  accurate  prediction  of  the 

steady  flow  field  is  essential.  However,  if  the  flow  field 
is  everywhere  subsonic,  i.c.,  <  (l-H«)/(  Y-t-DHi,  this 

dependency  is  not  as_iaportaat.vaincs..there  la  no  change  of 
equation  type.  -  -  - 


To  obtain  an  accurate  steady  flow  field,  particularly  in 
the  transonic  regime,  the  small  disturbance  theory  often  is 
inadequate  due  to  its  local  failure  at  the  leading  edge  and, 
in  some  cases,  to  the  limitation  on  shock  strength.  Although 
this  problem  can  be  circumvented  by  using  numerical  codes 
based  on  the  full  potential  equations,  or  even  better,  the 
Euler  equations,  the  cost  of  using  such  codes  in  engineering 
applications,  such  as  the  prediction  of  flutter  boundaries,  is 
prohibitively  high.  The  grid  systems  that  are  essential  for 
such  calculations  are  unnecessary  for  the  unsteady  flow  and 
may  even  be  harmful  due  to  stringent  time  step  requirements 
associated  with  locally  refined  meshes.  More  importantly, 
viscous  effects  have  been  ignored.  Following  Fung  and  Chung 
(4),  we  develop  a  method  that  allows  the  algorithm  developed 
here  to  be  applied  to  an  airfoil  for  which  we  have  an 
accurate  steady-state  pressure  distribution.  Thus,  this 
algorithm  can  be  Implemented  to  determine  the  unsteady 
response  of  an  airfoil  whose  steady  pressure  distribution  is 
known.  In  this  way,  at  least  steady  viscous  effects  are 
accounted  for.  As  we  shall  discuss  later,  it  is  much  easier 
to  get  accurate  steady-state  transonic  results  than  unsteady 
ones  because  experimental  flows  are  easily  contaminated  by 
the  reflections  of  unsteady  disturbances  from  wind  tunnel 
walls. 

Since  an  accurate  steady-state  pressure  field  is 
desired,  we  find  the  small  perturbation  flow  field  that 
corresponds  to  a  pressure  distribution  provided  from  wind 
tunnel  measurements  or  accurate  numerical  computations. 
This  is  referred  to  as  an  inverse  design  problem:  find  an 
airfoil  that  has  a  given  pressure  distribution. 

In  the  Inverse  problem,  since  pressure  is  given,  the 
values  of  on  the  slit  representing  the  airfoil  are  known, 
up  to  an  arbitrary  constant  c',  from  Integrating  To 

determine  this  constant  c',  a  closure  condition  is  required. 
In  small  disturbance  theory,  since  1>y  is  equivalent  to  the 
body  slope,  a  closed  airfoil  requires  that 

1 

/  ['^9(x,0'*')  -  'l'OCx,0~)ldx  -  0. 

0  ' 

Although  the  numerical  computations  cannot  normally  satisfy 
this  condition  exactly,  it  provides  a  means  of  determining 
the  constant  c'. 

An  Iterative  process  is  used  to  solve  Equation  (1)  for 
the  steady  state  The  value  Ac' found  from 

A  •  ^ 

If  -  «  {  /  [♦^(x.O'^  -  ♦f(x,0-)ldx-<}  (17) 

^  0 

Is  added  to  every  unsteady  value’  the  ntl^  iterate  of  <p  at 


until  a  steady  state  is  reached.  The  relaxation  factor  in 
Equation  (18)  helps  converge  Che  calculations  to  their  ■iC'^ady 
state;  the  small  parameter  c  is  added  to  yield  a  smooch  body 
slope  distribution.  Detailed  descriptions  of  this  procedure 
can  be  found  in  Fung  and  Chung  (4). 

Once  a  steady  flow  field  is  obtained,  the  unsteady 
response  is  calculated  using  the  algorithm  described  earlier 
or  Che  LTRAN2  code  developed  in  Ref.  |61  with  the  body  slopes 
determined  by  the  inverse  procedure. 

5.  INDICIAL  RESPONSE  METHOD 

A  major  advantage  of  time  linearization  is  the 
possibility  of  using  superposition;  chat  is,  adding  solutions 
chat  are  independent  to  form  a  composite  solution  in  time. 
Basically,  if  Che  system  used  is  linear  in  time,  including  the 
governing  equations  and  Che  associated  boundary  conditions, 
and  if  Pi(t)  and  solutions,  then  the  composite 

+  a2'^2  is  also  .a  solution.  For  example,  if  the  change  in 
lift  coefficient  of  an  airfoil  as  a  function  of  time  for  a 
step  change  in  angle  of  attack  is  Cj^^(t),  then  the  lift 
coefficient  ,Cj^(t),  for  a  variation,  a(c),  in  angle  of  attack  is 

Cj(t)  -  Cj(t)a(o)  -  J  C^^,(c’)  dt',  (18) 

o 

a  direct  application  of  Che  Duhamel  integral.  Through 
linearization,  the  system  of  Equations  (5),  (6a),  and  (6b)  is 
linear  and  its  dependence  on  k  can  be  removed  by  simply 
scaling  k.  Consequently,  sinusoidal  variations  of  any 
reduced  frequency  k  are  obtainable  from  the  indicia  1 
response  of  a  step  change  using  Equation  (18). 

6.  NUMERICAL  EXAMPLES 

To  demonstrate  the  merits  of  the  time  linearized 
method,  we  choose  a  set  of  measurements  for  a  NACA  64A010 
airfoil  at  Moo“  0.8  from  the  well  documented  experiments  of 
Davis  and  Malcolm  [7].  In  these  experiments,  Che  airfoil  was 
pitching  at  an  angle  l.D”  sinusoidally  over  a  range  of 
reduced  frequencies.  Figure  3  shows  the  experimental  steady 
pressure  distribution  on  this  airfoil.  This  pressure 
distribution  was  used  as  input  to  the  inverse  algorithm  IAF2 
[4]  that  provides  the  steady  flow  field  used  later  in  UTFC, 
Che  algorithm  based  on  Equation  (16).  Because  there  were 
nineteen  measured  pressure  values  from  these  experiments,  a 
100  X  80  grid  with  23  mesh  points  on  the  airfoil  was  chosen. 
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Figure  3.  Upper  and  lower  steady  surface  pressure  froa 
the  experiaents  of  Ref.  7  for  a  NACA  64A010  airfoil, 

-  0.8. 


Figure  4  compares  the  lift  and  moment  coefficients 
obtained  from  UTFC  with  those  from  the  experiments  and  from 
LTRAN2-HI,  a  nonlinear  code  based  on  Equation  f/a).  Overall 
agreement  between  the  measured  and  computed  values  exists. 
Note  also  how  well  the  results  from  UTFC  agree  with  chose 
from  LTRAN2-H1,  verifying  that  time  linearization  provides  an 
adequate  approximation  for  sufficiently  small  c. 
Furthermore,  it  seems  that  the  magnitudes  of  the  moment 
coefficient  from  UTFC  are  closer  to  Che  experimental  values 
chan  those  from  LTRAN2,  wtiich  is  attributed  to  the  better 
treatment  of  moving  shock  waves  in  the  time-linearized 
theory  when  the  amplitude  of  motion  is  small. 

Compared  to  Che  magnitudes,  which  are  in  excellent 
agreement  with  the  measured  values,  Che  phase  angles  are  in 
only  moderately  good  agreement.  Numerical  experiments 

discussed  in  the  next  section  indicate  chat  this  discrepancy 
is  due  to  wind  tunnel  Interference. 

It  is  worth  mentioning  that  ail  the  results  from  UTFC 
were  obtained  using  120  time  steps  per  cycle,  or  3“  of 
circular  oscillation  per  time  step,  for  all  frequencies  and 
that  no  instability  due  to  time  step  limitation  was  observed. 
However,  it  was  reported  by  Hessenlus  and  Goorjian  [8]  that 
in  order  to  produce  a  stable  solution  at  k  >  O.S  using 
LTRAN2-HI,  as  many  as  1440  time  steps  per  cycle  were 
necessary.  As  discussed  in  (9|,  this  Is  attributed  to  the 
leading  edge  over-expansion,  to  unnecessary  local  mesh 
refinements,  and/or  to  the  shock  over-shoot  as  a  result  of 
the  Murman-Cole  switching  that  may  admit  temporal  expansion 
shocks;  however  Goorjian  and  van  Busklck.  [9i,^also  reported 
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^^*®*^*  4.  Lift  and  leading  edge  aonent  coefficients  vs 
reduced  frequency  for  a  pitching  NACA  64A010  airfoil. 
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chat  the  monotone  scheme  of  Engquist-Osher  removes  the 
nonlinear  instability  of  over-expansion.  This  problem  does 
not  occur  in  the  time-linearized  formulation  because  the 
local  Mach  number,  which  determines  the  type  of  dif ferencin*; 
in  the  Murman-Cole  formulation,  depends  only  on  the  steady- 
state  values  and  this  is  unaffected  by  the  unsteady 

perturbations. 

7.  WIND  TUNNEL  WALL  EFFECTS 

The  accuracy  and  reliability  of  wind  tunnel  experiments 
have  been  difficult  to  assess,  ^particular ly  in  the  transonic 
regime.  The  presence  of  wind  tunnel  wails,  which  in  practice 
are  seldom  farther  than  ten  chords  away  from  the  testing 
airfoil,  makes  the  simulation  of  free  flight  conditions 
difficult.  Many  attempts  have  been  made  to  minimize  wall 
interferences,  including  wall  ventillation  through  slots  as 
used  in  the  Davis  and  Malcolm  |7J  experiments,  the  smart 
wind  tunnel  with  adjustable  wall  porosity  proposed  by  Sears 
(lOj,  and  the  flexible  wall  method  developed  by  Goodyer  (llj. 
However,  none  of  these  methods  were  designed  to  eliminate  or 
minimize  unsteady  wall  interference. 

Except  for  solid  walls,  where  the  boundary  conditions 
are  well  defined,  the  general  wall  conditions  in  experiments 
do  not  lend  themselves  to  analytical  modeling,  making 
quantitative  comparisons  between  numerical  and  experimental 
results  difficult,  if  not  impossible.  Figure  5  shows  typical 
variations  of  the  lift  coefficient,  C^,  as  a  function  of  time, 
measured  in  chord  length  travelled,  for  an  airfoil  after  a 
sudden  change  in  incidence.  Notice  that  after  a  lapse  of 
about  100  chord  lengths,  the  lift  coefficient  reaches  90%  of 
its  steady  value  -  the  solid  curve  computed  using  a  grid 
whose  boundaries  are  87  chords  away  from  the  airfoil  - 
indicating  that  the  wind  tunnel  wall  reflections  did  not 
have  enough  time  to  return  to  the  airfoil.  The  values  shown 
by  the  asterisks,  obtained  with  boundary  conditions 
satisfying  Equation  (4)  posed  at  the  grid  boundary  13  chords 
away,  agree  very  well  with  the  solid  curve,  implying  that 
Equation  (4)  is  a  good  representation  of  the  solution  at 
that  distance  and  that  no  significant  wall  interference  is 
detected  except,  perhaps,  some  minor  nonlinear  effects. 

Also  shown  in  Figure  3  are  the  results  (stars)  obtained 
using  the  non-ref lectlve  boundary  conditions  of  Engquist  and 
Hajda  [12]  and  Chat  using  the  constant  pressure  condition, 

♦x  ■  0* 

The  non-ref lective  boundary  condition. 


> 


stating  that  far  away  only  one  family  of  oucgola>;  w.ivc*s 
along  the  characteristics  dy/dx  »  achieves 

about  yOX  of  the  steady-state  lift  and  gives  a  subscantiil 
improvement  over  the  constant  pressure  condition,  which 
causes  a  30Z  reduction  in  magnitude  of  the  lift  coefficient 
and  a  shortened  characteristic  time  of  only  40  chords. 


Figure  5,  Indlcial  pitch  response  for  an  NACA  64A006  at 
**  0«88« 


It  is  observed  that  in  these  numerical  experiments  the 
lift  coefficient  reaches  a  new  steady  value  that  depends  on 
the  condition  at  the  far-fleld  boundaries  through  a  process 
describable  by  the  simple  model 

Ci(t)  ■  (20) 

where  1/X  denotes  a  characteristic  time.  Substituting 
Equation  (20)  into  Equation  (18)  and  assuming  harmonic 
pitching  motions,  a(t)  sin  kt,  the  unsteady  lift  coefficient 
is  then 

C£(t,  k)  - - -jjj  3in(kt  -  tan'^-j)  (22) 

*  <T>'j 

This  result  shows  chat  for  low  reduced  frequencies,  Che 
magnitude  and  Che  phase  depend  on  the  parameter  (k/X).  For 
t3rpical  transonic  flows  over  an  airfoil,  X  is  small  and 
decreases  rapidly  with  Mach  number.  Due  to  wall  confinement 
in  a  typical  wind  tunnel  experiment,  the  steady  value  C£(<») 
would  be  smaller  and  the  characteristic  time  l/X  would  also 
be  smaller  chan  those  measured  in  unbounded  flow.  Hence  the 
magnitude  of  the  unsteady  lift  coefficient  measured  in  a 
wind  tunnel  would  be  smalleY;;^'tCs  decrease  with  k  not  as 
rapid  due  to  the  larger  X.  Also'  the  phase  lag  increase  with 


k  would  not  be  as  rapid.  These  observations  are  consistent 
with  the  results  shown  in  Figure  k.  For  high  reduced 
frequencies,  the  effects  of  wall  boundary  condition  are 
poorly  understood.  In  the  case  of  a  solid  wall,  resonance  is 
known  to  occur,  but  most  tunnel  walls  are  ventilated, 
aiitlgating  the  effects  of  resonance. 

Figure  6  compares  the  same  set  of  experimental  values 
in  Figure  A  with  that  of  a  time  linearized  calculation  using 
the  steady  pressure  as  measured  in  the  experiment  as  input 
and  with  the  solid  wall  condition,  <{>y  =  0,  prescribed  at  the 
grid  boundary  7.2  chords  away  from  the  airfoil.  At  the 
frequencies  k  »  0.125  and  k  ■  0.325,  the  resonance  condition 
as  indicated  by  suddent  drops  of  magnitude  and  phase  seems 
to  occur.  Figure  7  shows  similar  comparisons,  but  with  the 
non-reflective  boundary  condition  described  earlier  applied 
at  7.2  chords.  For  moderate  reduced  frequencies,  k  >  0.2,  the 
computed  values  (triangles)  for  a  finite  boundary  as  close 
as  7.2  chords  away  are  the  same  as  those  for  an  unbounded 
region,  but  for  low  reduced  frequncies,  these  values  are 
closer  to  the  experimental  ones,  indicating  the  presence  of 
wall  interference.  Figure  8  shows  again  similar  comparisons 
with  the  linear  asymptotic  far-field  condition.  Equation  (4) 
prescribed  at  Xhe  same  distance,  but  with  the  phase  of  r(t) 
arbitrarily  adjusted  to  compensate  for  the  time  that  might 
have  taken  the  signal  to  reach  the  boundary.  Surprisingly, 
these  values,  especially  the  phase  angles,  coincide  with  the 
experimental  results,  Indicating  the  importance  of  wall 
interference  and  that  part  of  the  interference  is  due  to 
nonlinear  effects  not  accounted  for  in  a  strict  application 
of  Equation  (4). 

Although  a  general  description  of  wind  tunnel  wall 
interference  is  still  lacking,  these  numerical  experiments 
have  indeed  demonstrated  the  seriousness  of  the  problem. 
More  detailed  discussion  of  wall  interference,  including 
three-dimensional  effects,  can  be  found  in  Przybytkowski 
[131. 
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Flgvrc  7.  Coaparlsoa  of  the  lift  coefficients  using  the  non- 
refleetlee  boondsty  condition  st  s  distance,  j  •  7.2. 


CONCLUSION 


We  have  shown  in  this  chapter  that  f'jt  small  unsteady 
disturbances,  Che  time-linearization  technique  provides 
accurate  descriptions  of  the  flow  fields  despite  its 
simplicity.  It  allows  the  decoup  lint'  of  t!ie  unsteady  flow 
field  from  Che  steady  flow  field,  which  could  be  obtained 
more  efficiently  and  accurately  by  ocl\er  means,  includim’, 
direct  experimental  measurements.  The  extension  of  this 
procedure  to  three-dimensional  flow  fields  usin^  strip 
theory  for  high  aspect  ratio  wings  is  st  r  ii  ?lit  f  orw.ir  li,  th.- 
savings  are  substantial  and  the  success  warranted. 

The  numerical  examples  discussed  here  demonstrate  the 
importance  of  wall  interference  in  unsteady  transonic  winn 
tunnel  measurements,  which  deserves  a  better  understanding. 

ACKNOWLEDGMENTS 

This  work  was  supported  by  the  AFOSR  through  Grant  No. 
896875  and  by  NASA  through  Grant  No.  NGT-U3-002-80iJ. 

REFERENCES 

1.  LIN,  C.  C.,  T.  REISSNER,  H.  S.  TSIEN  -  On  Two-Dimensional 
Non-Steady  Motion  of  a  Slender  Body  in  a  Compressible 
Fluid.  J.  Math  Phys.,  Vol.  3,  pp.220-3l,  1948. 

2.  FUNG,  K.-Y.  -  Far  Field  Boundary  Conditions  for  Unsteady 
Transonic  Flows.  AIAA  Journal,  Vol.  19,  pp.  180-33,  1981. 

3.  YU,  N.  J.,  SEEBASS,  A.  R.,  BALLHAUS,  W.  F.  -  Implicit  Shock- 
Fitting  Scheme  for  Unsteady  Transonic  Flow  Computations. 
AIAA  Journal,  Vol.  16,  pp.  673-78.  1978. 

4.  FUNG,  K.-Y.  and  CHUNG,  A.  -  Computation  of  Unsteady 

Transonic  Aerodynamic  Responses  Using  a  Prescribed  Input 
Steady  State  Pressure  Distribution,  AIAA-82-0956. 

AIAA/ASME  3rd  Joint  Thermophysics,  Fluids,  Plasma  and 
Heat  Transfer  Conference,  St.  Louis ,  Missouri,  1982,  to 
appear  in  AIAA  Journal. 

5.  BALLHAUS,  W.  F.,  STEGER,  J.  L.  -  Implicit  Approximate- 

Factorization  of  Unsteady  Transonic  Flows  by  the  Indicia! 
Method.  AIAA  Journal.  Vol  16,  pp.l  17-24,  1975. 

6.  BALLHAUS,  W.  F.,  GOORJIAN,  P.  M.  -  Implicit  Finite- 

Difference  Computations  of  Unsteady  Transonic  Flows 
About  Airfoils.  AIAA  Journal,  Vol.  15,  pp.l728-35,  1977. 

7.  DAVIS,  S.  S.  and  MALCOLM,  G.  N.,  Experimental  Unsteady 

Aero-Dynamics  of  Conventional  and  Supercritical  Airfoils. 
NASA  TM-81221,  1980.  .  _ 


,  » 


*• 


8. 


HESSENIUS,  K.  A.  and  GOORJIAN,  P.  M.  -  Validation  of 
LTRAN2-H1  by  Comparison  with  Unsteady  Transonic 
Experiment.  AL\A  Journal,  Vol.  20,  pp.  731-32,  ann  NASA 
TM-81037,  1982. 


9.  GOORJLAN,  P.  M.  .and  Van  Buskirk,  R.  -  Impli.cic  Calculations 
of  Transonic  Flows  Using  Monotone  Methocls.  A  LA  A  J9ti', 

Aerospace  Sciences  Meeting  Sc.  Louis,  Missouri,  19ai. 


lO.  SEARS,  W.  R.  -  A  Note  on  Adaptive-Wail  Wind  T.in.ielj, 
Journal  of  Applied  Mathematics  and  Physics,  Vol.  J.h,  pp. 
91  5-27,  1977. 


11.  GOODYER,  M.  J.  -  The  Self-Streaming  Wind  Tunnel,  NASA 
TMX-72699,  1975. 

12.  ENGQUIST,  B.  and  MAJDA,  A.  -  Numerical  Radiation  Boundary 
Conditions  for  Unsteady  Transonic  Flows.  Journal  of  Com¬ 
putational  Phvslcs,  Vol-  40,  pp,  91-103,  1981. 

13.  PRZYBYTKOWSKI,  S.  M.  -  Effect  of  Wall  Interference  in 
Unsteady  Transonic  Flows.  Ph.D.  Thesis,  University  of 
Arizona,  1983. 


Vorticity  at  the  Shock  Foot  in  Inviscid  Flow 

K.-Y.  Fung 


VokMiw  21 ,  Numbr  6.  Jum  1983.  Pag*  91S 


■  'fi  H.  nmer  Technical  Librar, 

J1333 

!’"-=irr,h  I.3hor’'--- 


AMEMCAN  INSTITUTE  OF  AERONAUTICS  ANO  ASTRONAUTICS  *1290  AVENUE  OF  THE  AMERICAS  •NEW  YORK,  NEW  YORK,  N.Y.  10104 


JUNE  1983 


TECHNICAL  NOTES 


915 


proximation  is  made,  and  this  is  inconsistent  with  the  con¬ 
servation  of  normal  momentum  across  the  shock.  Also,  only 
the  weak  form  of  the  governing  equation  is  satisfted  to  an 
order  normally  so  restricted  by  numerical  stability  that  the 
actual  shock  is  smeared  over  several  grid  points.  As  the  shock 
gets  stronger,  entropy  and  vorticity  production  behind  the 
shock  can  no  longer  be  ignored  and  potential  theory  fails. 
These  effects,  added  to  the  already  complex  flowfield,  makes 
the  construction  of  a  proper  numerical  scheme  a  difficult 
task. 

In  this  Note  we  study  the  flowfield  immediately  down¬ 
stream  of  a  shock  at  its  root  where  the  shock  meets  a  smooth 
convex  surface.  Lin  and  Rubinov'  Tirst  noted  that  a 
singularity  occurs  at  the  shock  root.  They  also  argued  that  the 
shock  shape  at  the  root  must  be  of  the  form 

f =V'^ 
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IntrtMiuction 

IT  is  characteristic  of  transonic  flows  to  have  a  shock  or 
shocks  embedded  in  the  flowHeld.  The  flow  immediately 
behind  the  shock  is  related  to  the  flow  ahead  of  it  by  the 
Rankine-Hugoniot  conditions  and  is  a  function  of  the  shock 
shape.  If  the  shock  is  normal  to  the  body  surface,  the  flow 
behind  the  shock  will  be  subsonic  and  its  shape  will,  in 
general,  not  be  known  a  priori.  The  shock  shape  must  be 
determined  in  conjunction  with  the  local  flowfield. 

Although  inviscid  transonic  flow  past  a  body  can  be 
computed  routinely  using  numerical  methods  developed  in  the 
last  decade,  the  shock  region  has  always  been  the  most 
erroneous  part  of  the  solution.  Often  the  potential  ap- 
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where  (,  n  are  coordinates  of  the  shock  measured  along  and 
normal  to  the  body,  respectively.  Zierep^  also  found  the  same 
shock  shape  but  was  unable  to  determine  the  constant  k  for  a 
convex  body.  Later  Gadd’  pointed  out  that  the  flow  behind 
and  at  the  shock  root,  determined  by  the  Rankine-Hugoniot 
conditions,  experiences  a  discontinuity  in  curvature  in  order 
to  conform  to  the  body.  Such  a  flow  is  known  to  have  a 
multivalued  normal  pressure  gradient  and  a  streamwise 
pressure  gradient  that  is  logarithmically  singular. 

We  shall  determine  the  vorticity  behind  the  curved  shock  at 
this  singular  point  and  discuss  to  what  extent  the  flowfield  is 
affected  by  this  vorticity. 

Shock-Induced  Vorticity 

It  is  well  known  that  the  vorticity  behind  a  shock  can  be 
computed  by  applying  Crocco’s  theorem,  i.e., 

f  -  ‘  f  J  dgj  /  dP;.  ] 

<7,sin((r-a)  12  df  p,  df  J  ' 

where  subscript  2  denotes  quantities  evaluated  after  the 
shock,  t  is  the  vorticity  induced  by  the  shock,  q  the  flow 
speed,  p  the  density,  p  the  pressure,  f  the  distance  along  the 
shock,  a  the  shock  angle  measured  relative  to  the  upstream 
flow,  and  a  the  flow  deflection  angle  after  the  shock  (Fig.  I). 

The  Rankine-Hugoniot  conditions  require  that  the  post¬ 
shock  quantities  be  related  to  preshock  quantities  as  follows: 


Pi=P;+P/<7)(/-Osin'o 

Py  =  IZi  + 

P2  I'ki  (i  + 


a 
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where  the  subscript  /  denotes  quantities  ahead  of  the  shock,  y 
is  the  ratio  of  the  specific  heats,  and  M  is  the  Mach  number. 
After  substituting  Eq.  (2)  into  Eq.  (I),  we  see  that  the  vorticity 
equation  becomes 


rsino 


1 


(f  (A/J  -  /)sin^o-cosal 


/ 

/>;<?( 


dP/ 

df 


-  (/-f)sinacosa 


do ' 
Jtl 


(3) 


At  the  shock  root  where  the  shock  is  normal  this  formula  can 
be  further  simplified  to 


(f-Oq, 

f  L  R 


(J-e)cosa 


da 

df 


] 


(4) 


where  R  is  the  radius  of  curvature  of  the  body  at  the  shock 
root.  We  note  that  the  last  term  on  the  right-hand  side  of  Eq. 
(4)  cannot  be  evaluated  immediately  but  cosuldu/df)  must  be 
negative  since  sino  attains  its  maximum  at  the  root  with 
o  =  t/2.  We  thus  conclude  that  the  vorticity  is  clockwise  or 
negative,  and  is  a  second-order  quantity  proportional  to  (I  - 
e)^  since,  from  Eq.  (2), 

/  _ ,  =  2{M]sin^a-l) 

(■y-k  /lAfJsin-'o 

In  order  to  evaluate  the  vorticity,  the  shock  shape  a(f)  must 
be  found.  We  note  here  that  the  potential  approximation  is 
only  valid  to  second  order  at  the  shock  root  unless 
1/R  =  0(I  -(),  i.e.,  unless  the  slender  body  approximation  is 
made. 


can  be  very  informative  in  understanding  the  flow  behind  a 
shock. 


Vorticity  at  the  Shock  Root 

Since  vorticity  is  shown  to  be  of  the  order  of  (I  a 

formal  expansion,  e.g.,  ♦  =  V“-k5V' -k...  in  a  small 
parameter  6=6(1  -f ),  yields  the  lowest  order  equation  for  ♦ 
valid  at  the  vicinity  of  the  shock  root  as  follows 


(9) 


The  solution  of  this  equation  subject  to  the  stated  boundary 
conditions  was  first  given  by  Gadd^  and  later  solved  in  a  more 
formal  manner  by  Oswalitsch  and  Zierep.’  The  shock  shape 
they  found,  which  is  different  from  that  predicted  by  Lin  and 
Rubinov,'  is  normal  to  the  body  and  has  a  logarithmic 
singularity  in  curvature  at  the  body.  The  unknown  quantity, 
coso(do/df),  evaluated  at  the  shock  root  is  then  zero. 

We  conclude  that  since  vorticity  is  second  order  the 
flowfield  behind  the  she  ck  remains  irrotational  to  the  lowest 
order.  The  local  curvature  of  shock  at  the  shock  root,  despite 
being  logarithmical-infinite,  does  not  contribute  to  the 
vorticity.  For  slender  bodies  the  shock  curvature  as  well  as  the 
body  curvature  are  0(1 -r).  Thus  our  result  for  vorticity 
reduces  to  the  usual  small  disturbance  result  of  (I  -t)^. 

Since  all  the  higher  order  equations  are  of  Poisson-type 
with  homogeneous  boundary  conditions,  their  solutions 
should  be  regular  and  have  no  more  singular  contribution  to 
the  shock  curvature  at  the  root  than  the  lowest  order 
logarithmic  one.  We  then  obtain  the  equation  for  vorticity  at 
the  root  as 


Flow  Behind  (he  Shock 

The  inviscid  subsonic  flow  after  the  shock  would  be 
rotational  in  general  and  therefore  is  governed  by  the  Euler 
equations.  The  continuity  equation,  written  in  intrinsic 
coordinates  s  and  n,  along  and  normal  to  flow,  is 


h=--'  (10) 

where 

e=  -k  ^ 
y+l  (y  +  J)Mj 


9pg 

9s 


(6)  in  terms  of  known  quantities  q,  and  M,,  upstream  of  the 
shock  and  the  body  curvature  \/R. 


The  vorticity  equation  by  definition  is 

^  39  9q 
^  9s  9n 


P) 


where  9  is  the  flow  angle  measured  with  a  fixed  reference 
frame. 

If  one  defines  an  intrinsic  stream  function  V  as 


Concinsion 

We  have  shown  that  vorticity  behind  the  shock  is  of  the 
order  of  (I  -  f )  ■^  even  at  the  shock  root  where  the  curvature  of 
the  shock  is  infinite,  and  have  obtained  a  formula  for  this 
vorticity  valid  at  the  root  in  terms  of  upstream  quantities 
only.  This  formula  can  be  used  in  the  construction  of  an 
accurate  numerical  scheme  for  Eq.  (8). 


4',  =f<iP9  and  ♦,  =  -9 

then,  Eq.  (6)  is  automatically  satisfied  and  Eq.  (7)  becomes 
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ll+(y-I}Afn  ^ 
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(8) 


This  equation  can  be  solved  for  with  given  fas  in  Eq.  (4)  by 
specifying  tun^0;t,  and  4'„  at  the  shock  boundary  f(s,n) 
with  4',  and  4',  being  finite  at  downstream  infinity.  The 
requirement  that  this  elliptic  system  is  not  being  over-specified 
thus  determines  the  shock. 

In  general,  Eq.  (8)  can  only  be  solved  by  numerical  means; 
however,  a  local  expansion,  making  use  of  the  fact  that  the 
vorticity  is  only  a  second-order  quantity  as  mentioned  earlier. 
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